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AN  APL  FUNCTION  FOR  BIVARIATE  NORMAL  PROBABILITIES 


INTRODUCTION:  Bivariate  normal  distributions  have  many 
applications  such  as  in  combat  modelling,  weapons  systems 
effectiveness  studies  and  weather  prediction  problems;  several 
other  applications  are  discussed  in  [4].  It  is  well  known  that 
probability  statements  for  a  general  bivariate  normal 
distribution  can  be  transformed  into  equivalent  statements  for  a 
standard  bivariate  normal  distribution  with  0  means  and  variances 
equal  to  1.  It  is  thus  sufficient  to  be  able  to  compute 
cumulative  probabilities  for  a  standard  bivariate  normal 
distribution  with  probability  density  function 


1 


2(l-p2) 


|_x  -  2pxy  +  y 


f  (X,y)         =  /  I  e  _00<X<0O_0O<y<0O 

2  TTV 1  -  P 

Similar  to  the  case  with  the  univariate  normal  distribution,   the 
bivariate  cumulative  distribution  function  (c.d.f.) 


k  h 
F(h,k)  P[X<h,Y<k]         f   f   f(x,y)   dx  dy 


I  I 

-00  —00 


does  not  have  a  closed  form  solution  and  numerical  integration  is 
the  usual  approach  to  evaluating  such  integrals.  Because  of  the 
importance  of  the  distribution,  the  National  Bureau  of  Standards 


(NBS)  has  published  an  extensive  set  of  tables  for  P[X>h,Y>k] 
[5].  Another  set  of  tables  using  a  different  approximation,  has 
been  generated  by  Owen  [6].  Other  ways  to  approximate  the 
bivariate  normal  integral  are  discussed  in  [1],  [3]  and  [7].  All 
of  these  approximations  involve  numerical  integration  and  require 
extensive  computer  programming,  rendering  them  to  be  not  readily 
suitable  for  obtaining  on  the  spot  results.  These  days,  with  the 
easy  availability  of  microcomputers,  it  would  be  useful  to  have  a 
program  to  compute  bivariate  normal  probabilities  interactively 
and  also  be  able  to  incorporate  such  a  program  within  other 
application  programs.  This  will  allow  an  analyst  to  perform 
sensitivity  studies  by  varying  the  input  parameters  in  an 
application  program  and  observing  the  effect  on  the  measures  of 
effectiveness  of  interest. 

Recently,  Wang  [8]  developed  a  new  algorithm  to  compute 
bivariate  normal  probabilities,  that  does  not  require  numerical 
integration,  relatively  easy  to  program  and  provides  quite 
accurate  results.  Investigations  by  Wang  indicate  that  the 
computed  probabilities  compare  very  well  with  those  in  the  NBS 
tables  and  the  computer  resources  needed  are  not  excessive.  The 
approach  used  by  Wang  is  to  start  with  an  approximate  contingency 
table  of  cell  probabilities  for  a  rectangular  grid  on  the  x,  y 
plane  and  then  apply  the  "iterative  proportional  fitting 
algorithm"  [3;  sec  3.5]  to  modify  the  cell  probabilities;  the 
iteration  is  continued  until  the  marginal  probabilities  in  the 
contingency  table  coincide  with  their  true  values  to  within  a 
specified  degree  of  accuracy.  Since  the  marginal  distributions 
of  the  random  variables  X  and  Y  are  univariate  normal,  the   exact 


marginal  probabilities  can  be  determined  using  computer  programs 
available  in  most  statistical  software  packages  or  by  writing  a 
subprogram  for  the  purpose.  It  should  be  noted  that,  although  a 
bivariate  normal  distribution  is  defined  over  the  entire  plane, 
for  all  practical  purposes  it  is  sufficient  to  consider  only  the 
square  subregion  [-4,4]  x  [-4,4]  since  the  probability  content 
outside  of  this  subregion  is  negligible. 

Wang's  algorithm  is  defined  by  the  following  steps: 

1.  Let  -4  =  aQ   a^    ...    am  =  4  be  a  partition  of 
the  interval  [-4,4]  along  the  x  -  axis  and  -4  =  bQ  b-^ 

.  .  .    bn  =  4,  a  partition  along  the  y  -  axis;  identical 
equal   length  partitions   for  both  axes   is   recommended  to 
reduce  computational  complexity. 

2.  Let      Ax  =  Xj_  -  aj.-l  i  =  ■*-'  2,...,  m 

Ay  =  bj  -  bj_x    j  =  1,  2,...,  n 
m^  =  a:  i  +  a^    i  =  l,2...,m 

n^  =  b.:  n  +  b.;     j  =  l,2,...,n 


p  =   p    (1  -  Ax^)   (1  -  Ay2  ) 

12         12 
where    is  the  correlation  coefficient  (specified) 

3.   Let   p.  =  p  [a.  .  <  x  <  a.]   L  =  1,2,..., in 
i       i-l     —  l 

p.  =  P  [b.    <  Y  <  b.]        j  =  1,2,. ...n 
J       j-l     -  J 


be  the  marginal  probability  contents  of  the  i-th  and  j -th 
subintervals  along  the  x  and  y  axes  respectively. 


4.   Compute   p^j  ,  the  starting  approximate  probability  of  the 
(i,  j)  th  cell  as 


p  m.n. 
i  3 


^ 


=  e 


i-P2 


i  —  1/2,  ...m;  j  =  1,2,.  ..,11 


5.  The  application  of  the  iterative  proportional  fitting 
algorithm  results  in  the  following  equations  for  the  modified 
probability  of  the  (i,j)  th  cell  after  the  k-th  iteration: 
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k  =  1,3,5, 


k  =  2,4,6, .. . 


6.   Continue  the  iteration  process  until  for  some  even  number  k 


.00 


P. .   -  P. 
ij     i 
3=1 


<   £ 


and 


m 


i=l 


00   _  p 
ij 


<   £ 


where  £  is  a  prespecified  degree  of  accuracy  with  which  the  true 
marginal  probabilities  agree  with  the  marginals  in  the 
contingency  table. 


7.   To  compute   P  [  X  <  h  ,  Y  <  k  ]  (or  P  [  X  >  h,  Y  >  k  ]  )   sum 
the  probabilities  in  the  contingency  table  over  those  cells  for 

which  a  <  h  (a  >  h)  and  b .  <  k  (b .  >  k) .  In  those  cases  where  either 
i  —     i  J  —     3 

h  #  a.  and/or  k  ±    b.  for  any  of  the  partition  points  a±   and  b_.  ,  the 


accuracy  of  the  approximation  can  be  improved  by   including  h 
and/or  k  as  additional  partition  points. 

An  APL  function,  called  BVN,  to  compute  bivariate  normal 
probabilities  (both  P  [  X  <  h  ,  Y  <  k  ]  and  P  [  X  >  h,  Y  >  k  ])  is 
presented  in  the  appendix.  This  function  invokes  another  APL 
function  called  NCDF  to  compute  the  marginal  univariate  normal 
probabilities.  The  BVN  function  can  be  run  on  an  IBM-PC  /  AT 
compatible  microcomputer  using  an  APL  language  system  such  as 
APL*PLUS  from  the  STSC  corporation;  the  function  can  also  be  run 
under  VSAPL  on  the  NPS  mainframe  computer.  The  function  runs 
interactively  and  calls  for  keyboard  input  of  the  desired  equal 
length  partition  size  for  the  x  and  y  axes  and  the  degree  of 
accuracy  e  in  approximating  the  marginal  cell  probabilities. 
With  only  minor  modifications,  the  function  can  be  imbedded 
within  another  APL  function  as  a  subprogram.  Computations  using 
the  BVN  function  indicate  that  with  partition  size  x  =  y  =  .2 
for  both  the  x  and  y  axes  and  e  =  .00005,  a  four  decimal  place 
accuracy  as  compared  with  the  NBS  tables  can  be  achieved.  The 
computational  time,  as  is  to  be  expected,  increases  with  a 
decrease  in  the  partition  size  x  or  y,  an  increase  in  and 
to  a  lesser  degree  a  decrease  in  e  .  With  x  =  y  =  . 2 , 
e  =  .00005  and  .  1<  <.8  the  computational  time  (clock  time) 
was  between  90  and  180  seconds  on  the  Zenith  Z-248  (an  AT  type) 
microcomputer,  and  on  the  IBM  3033  mainframe  computer  these  times 
were  between  1  and  20  seconds.  The  computational  time  can  be 
reduced  considerably  (to  about  30  seconds  on  the  Z-248)  by 
choosing  x  =  y  =  .5  but  then  only  a  three  decimal  computational 


accuracy  can  be  expected.  For  a  fixed  value  of  if  the  c.d.f. 
is  to  be  calculated  for  several  choices  of  (11,10,  the  BVN 
function  needs  to  be  run  only  once;  with  a  very  minor 
modification  the  contingency  table  of  bivariate  cell 
probabilities  can  be  saved  in  a  matrix  and  all  that  is  left  is  to 
sum  the  probabilities  of  the  appropriate  cells.  If  needed,  it 
would  be  quite  straight  forward  to  generate  tables  for  various 
choices  of   and  (h,k) . 

Professor  I.  O'muircheartaigh  of  the  O.R.  department  and  I 
are  in  the  process  of  completing  a  Fortran  program  for  the 
problem  and  expect  to  submit  the  code  for  publication. 
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APPENDIX 


THIS  APPENDIX  CONTAINS  THE  LISTING  OF  TWO  APL  VARIABLES  BVNHOW  AND 
NCDFHOW  AND  TWO  APL  FUNCTIONS  BVN  AND  NCDF.  THE  TWO  HOW  VARIABLES 
PROVIDE  SHORT  DESCRIPTIONS  OF  THE  COMPUTATIONAL  SCHEMES,  THE  INPUT 
PARAMETERS  AND  THE  SYNTAX  FOR  THE  TWO  FUNCTIONS. 

BVNHOW 


THE  FUNCTION  BVN  COMPUTES  THE  C.D.F.  OF  A  STANDARD  BIVARIATE 
MORMAL  DISTRIBUTION  WITH  CORRELATION  COEFFICIENT  f> ,  USING  AN 
ALGORITHM  PROPOSED  BY  YUCHUNG  J.  WANG  (BIOMETRIKA,  19B7,  NO. 74,  185-90) 
THE  ALGORITHM  CONSISTS  OF  PARTITIONING  THE  X-Y  PLANE  INTO  RECTANGULAR 
:ELLS,  AN  INITIAL  APPROXIMATION  OF  THE  CELL  PROBABILITIES  AND  AN 
ITERATIVE  SCHEME  TO  MODIFY  THE  CELL  PROBABILITIES.   THE  ITERATION 
PROCESS  IS  TERMINATED  AS  SOON  AS  THE  MARGINAL  PROBABILITIES  COINCIDE 
giTH  THEIR  EXACT  VALUES  (THAT  ARE  UNIVARIATE  STANDARD  NORMAL 
PROBABILITIES, COMPUTABLE  USING  THE  APL  FUNCTION  NCDF)  TO  WITHIN  A 
SPECIFIED  DEGREE  OF  ACCURACY  e  .   THE  SYNTAX  FOR  THE  FUNCTION  IS 

RHO  BVN  W 

IHERE  RHO  IS  THE  CORRELATION  COEFFICIENT  AND  W  =  <x,y)  IS  THE  POINT 
\T   WHICH  THE  C.D.F.  IS  TO  BE  COMPUTED.   THE  FUNCTION  WILL  CALL  FOR 
HE  INPUT  OF  THE  DESIRED  LENGTH  FOR  AN  EQUAL  PARTITIONING  OF  THE 
NTERVAL  C*4,4:  ALONG  THE  X  AND  Y  AXES  AND    e    THE  DESIRED  DEGREE 
IF  ACCURACY  IN  THE  MARGINAL  PROBABILITIES.   THE  RECOMMENDED  CHOICES 
RE   .2   FOR  THE  PARTITION  LENGTH  Ax,  AND  .00005  FOR   s  .   THE  TOTAL 
OMPUTATION  TIME,  ON  THE  ZENITH  Z-248  MICROCOMPUTER,  SHOULD  BE  BETWEEN 
0  AND  ISO  SECONDS  DEPENDING  ON  THE  SIZE  OF  RHO. 

NCDFHOW 


HE  FUNCTION  NCDF  COMPUTES  THE  C.D.F.  OF  A  STANDARD  NORMAL  DISTRIBUTION, 
SING  THE  APPROXIMATION  DEFINED  IN  EQUATION  26.2.17,  PAGE  932  IN 
HE  HANDBOOK  OF  MATHEMATICAL  FUNCTIONS,  M.ABRAMOWITZ  AND  I.A.STEGUN 
DITORS,  PUBLISHED  BY  THE  NATIONAL  BUREAU  OF  STANDARDS  (REF.  C5]>.   THIS 
^PROXIMATION  IS  ACCURATE  TO  ATLEAST  TO  7  DECIMAL  PLACES.  THE  SYNTAX 
IDR  THE  FUNCTION  IS:    NCDF   Z      WHERE  Z  IS  AN  INCREASING  ARRAY 
lr  NUMBERS  FOR  WHICH  THE  C.D.F.  IS  TO  BE  COMPUTED.   THIS  FUNCTION  IS 
EVOKED  BY  THE  BVN  FUNCTION  TO  COMPUTE  MARGINAL  PROBABILITIES. 

-9- 


7  RHO  BVN  W;X;Y;RBAR;R;A;M;B;E;Q;D; I; J; K; LI ; L2; Al ; A2; Ml ;M2;P1 ; P2 
C13    '  ' 
C23    '  ' 

LZ1  ft THIS  FUNCTION  COMPUTES  THE  CUMULATIVE  PROBABILITIES   OF  A  STANDARD 

[43  fiBIVARIATE  NORMAL  DISTRIBUTION  (  means  0  and  st.devs  1)  WITH 

C53  ft CORRELATION  COEFFICIENT  RHO.   THE  SYNTAX  FOR  THE  FUNCTION  IS: 

C63  ft RHO  BVN  W   ,   WHERE   N-<x,y).   THE  FUNCTION  FIRST  COMPUTES  A 

C73  ft CONTINGENCY  TABLE  OF  CELL  PROBABILITIES  OVER  A  USER  DEFINED  PART IT 10 

C8D  ft OF  THE  X  AND  Y  AXES  AND  THEN  CUMULATES  THE  PROBAILITIES  IN  THE 

C9D  ft APPROPRIATE  CELLS.   THE  USER  IS  PROMPTED  TO  INPUT  THE  LENGTH  OF  THE 

C103  ftSUBINTERVALS  IN  THE  PARTITION  OF  ("4,4)   (e.g.,  .05)  AND  THE  DESIRED 

C113  ft  ACCURACY  <e.g.,  .0005  FOR  A  3-DECIMAL  ACCURACY).   THE  FUNCTION 

C123  flOUTPUTS  BOTH  Pr  C  X  <  x  ,  Y  <  y  3  AND  Pr  C  X  >  x  ,  Y  >  y  3. 
C133   '  INPUT  THE  DESIRED  PARTITION  SUBINTERVAL  LENGTH  FOR  X  AND  Y  AXES' 

C143  K<-0 
C153   '  ' 

C163   'INPUT  THE  DESIRED  ACCURACY  OF  COMPUTATIONS' 

C173  E«-0 

C1S3  X+WCU 

C193  Y+WC23 

C203  -►ENDlxv(XS-4)vYS-4 

C213  +END2>a(X>4)'NY>4 

C223  +END3XT.  (X>4)vy>4 

C233  X*L/4,X 

C243  Y*L/4,Y 

C  25  3  RBAR«-RHOx  1  -  ( K»2 )  -5- 1 2 

C  26 1  R*RBAR-s-  ( 1  -RBAR*2 ) 

C273  A<-"4+0,Kxt,S-5-K 

C283  AH-((X>A)/A)  fX,  (X<A)/A 

C293  A2*((Y>A)/A) ,Y, (Y<A)/A 

C303  Ml+(  (liAl)+-liAl)-r2 

C313  M2«-(  (1  +  A2)+-1  +  A2)*2 

C323  Ll+pMl 

C333  L2+PM2 

[34  3  B<-*RxMlo.xM2 

C353  P1«-(NCDF  1  +  AD-NCDF  "1  +  A1 

C363  P2«-(NCDF  UA2)-NCDF  *1  +  A2 

C373  REPEAT:  Q«-P1*+/B 

C383  D<-((Ll,Ll)|£>Q)x(  (LI  ,L1 )  pi  ,LlpO) 

C393  B+D+.xB 

C403  Q+P2++^B 
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[413   D<-<<L2,L2)|0Q)x<<L2,L2)  pl,L2pO) 
[423   B«-B+.xD 

[433   +REPEATXT.  ((+/<!  P1-- »-/B>  >E)  >0>  v  (+/  <  |  P2-+/B)  >E)  >0 
[443   I«-+/X>A1 
[453   J*+/Y>A2 

[463  'Pr  [  X  <  ',<*X>,'  ,  Y  <  ' , <5Y> , '  3  = 
[473  'Pr  [  X  >  ',(5X),'  ,  Y  >  '  ,  (*Y)  ,  '  3  - 
[483   +0 

[493  END1:  '  Pr  [  X  <  '  ,  <*X>  ,  '  ,  Y  <  '  , <?Y> , 
C503   '  Pr  [  X  >  ',(5X),'  ,  Y  >  ' , (*Y> , '  3 
[513   +0 

1523  END2: '  Pr  [  X  <  ',<5X),'  ,  Y  <  ',(?Y), 
:533   '  Pr  [  X  >  ',(5X),'  ,  Y  >  '  ,  (SY)  ,  '  3 
!543   +0 

1553  END3: '  Pr    [  X  <  ',<*X>,'  ,  Y  <  ',(5Y), 
!563   '  Pr    I    X  >  ',(5X),'  ,  Y  >  '  ,  <?Y)  ,  '  3 
573 
7 


'  ,3+/+/<I,J)tB 
' ,*+/+/<!, JJ+B 


1  ' 


0' 


1  ' 


«   0' 


',*<NCDF  X)xNCDF  Y 
',5<1-NCDF  X)x(l-NCDF  Y) 
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7  R+NCDF  Z;B;p;z;T;N;M;P 

[13  A  THIS  FUNCTION  COMPUTES  THE  C.D.F.  Pr  [  Z  <  z  ]  OF  A  STANDARD 

[23  ANORMAL  DISTRIBUTION  USING  THE  APPROXIMATION  IN  EQUATION  26.2-17, 

[33  A PAGE  932  IN  THE  HANDBOOK  OF  MATHEMATICAL  FUNCTIONS,  EDITED  BY 

C43  RABRAMOWITZ  AND  STATGUN,  NATIONAL  BUREAU  OF  STANDARDS. 

[53  Z<-,Z 

L6i  z+zc*za 

[73  N<-/»<Z<0)/Z 

cs:  m^(z>o)/z 

[93  Z*IZ 

[103  p*0. 2316419 

[113  T^-5-l+pxZ 

[123  z*<*-(Z*2)-i-2>-s-<o2)*0.5 

[133  B*  0.31 938 153  "0.356563782  1.781477937  "1.821255978  1.330274429 

[143  P^-1-2X(To.*t,5)+.xB 

[153  R+U-N+P)  ,  <-M)tP 
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